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Abstract 

We propose a hybrid model, coupling Lattice Boltzmann and Molecular Dynam- 
ics models, for the simulation of dense fluids. Time and length scales are decoupled 
by using an iterative Schwarz domain decomposition algorithm. The MD and LB 
formulations communicate via the exchange of velocities and velocity gradients at 
the interface. We validate the present LB-MD model in simulations of flows of liq- 
uid argon past and through a carbon nanotube. Comparisons with existing hybrid 
algorithms and with reference MD solutions demonstrate the validity of the present 
approach. 

1 Introduction 

The advent of nanotechnology provides us today with nanoscale devices that can affect 
flow phenomena that are important for several technological applications. Examples in- 
clude microffuidic channels with nanopatterned walls , biosensors embedded in aqueous 
environments [TJ |21 E3 E| , or bluff bodies with superhydorphobic surfaces 

The detailed investigation of flow physics at the nanoscale has been pioneered by Koplik 
[H] using Molecular Dynamics (MD) models. When nanoscale devices are active parts 
of micro and macroscale systems, a hybrid approach is required to integrate atomistic 
simulations with computational methods suitable for larger scales. A number of hybrid 
models coupling atomistic to continuum descriptions of dense fluids have been proposed |Hl 

EE EH unman!. 

O'Connell and Thompson [8 a coupled an atomistic with a continuum system where the 
average momentum of the overlap particles is adjusted through a boundary force. Flekk0y 
et al. [9. presented a model based on direct flux exchange between atomistic and continuum 
regions. Hadjiconstantinou and Patera |XD] proposed to decouple time scales by using the 
Schwarz domain decomposition method coupling an atomistic to a continuum description 
of a fluid. Convergence to a steady state solution is achieved through alternating iterations 
between steady state solutions within the atomistic and continuum subdomains. Nie et 
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al. |TT1 IT2] employed a domain decomposition algorithm to simulate Couette flow over 
a nanopatterned surface, as well as a cavity flow where the singularity at the corners 
between static and moving walls was described atomistically. In these simulations a singular 
boundary force is employed in order to compensate for the elimination of periodicity in 
the MD system. Werder et al. |T3] proposed an algorithm using the alternating Schwarz 
method to couple an MD model to an incompressible Navier-Stokes (NS) finite volume 
solver. They proposed a novel boundary force, based on the physical characteristics of 
the fluid that is being simulated. They reported on simulations of flows past a CNT and 
noticed average departures from reference MD simulations of the order of 4%. 

Here we extend the model proposed in by replacing the finite volume solver that 
was provided by a commercial software package (STAR-CD) by a Lattice Boltzmann (LB) 
method solving the incompressible NS equations. The proposed extension aims to 
take advantage of the mesoscopic modeling inherent in LB simulations and to allow for a 
broader geometric flexibility than the one allowed by the Finite Volume solver. In addition, 
in the present work we enhance the exchange between atomistic and continuum domains 
by not only communicating velocities as in ^3] but also velocity gradients. 

The paper is organized as follows. In section EJ we describe the hybrid model by 
outlining the MD and LB methods and describe their coupling. Results of liquid argon 
flows past and through a CNT are presented in section El We compare the flows to the 
reference MD solutions and discuss the computational efficiency of the hybrid model and 
conclude in section HI 



2 The LB-MD model 

In the present hybrid method an MD model describes the flow in the vicinity of a carbon 
nanotube (CNT) while a LB approach is used to simulate the behaviou of the continuum 
system away from the CNT. 

2.1 Molecular dynamics 



The atomistic region is described by MD simulations where the positions = (xi, yi, Zj) 
and velocities Uj = (u Xj i,u y) i,u Z) i) of the particles evolve according to Newton's equations 
of motion 



~dt 



at 



Evu(. 



(i) 



where F, and rrii are the force and mass of particle i, and is the distance between the 
particle and rj. Here we consider Lennard- Jones (LJ) model of argon interacting with 
CNTs. The potential Ufrij) is defined as 
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where €ab and gab are energy and length parameters, A and B denote species. The 
LJ interaction parameters for argon-argon and argon-carbon interactions are respectively 
^ArAr = 0.9960 kJmol" 1 , OArAr = 0.3405 nm, €ArC = 0.5697 kJmol" 1 , gatC — 0.3395 nm. 
The term Ub is the boundary potential and accounts for the interaction of the boundary 
region with the surrounding medium. It is further described in The CNT is modeled 
as a rigid structure to facilitate the investigation of the flow of argon. All interaction 
potentials are truncated for distances beyond a cutoff radius r c — 1.0 nm. The equations 
of motion (j2J) are integrated using a leap-frog scheme with a time step 5t = 10 fs. 

A desired velocity is enforced by relaxing the center of mass velocity Uk = l/Afc Y^iek u « 
of the Nk particles within a cell k towards according to a parameter A. The velocities 
Uj of the particles in the cell k are updated as 

U, = U; + X(u d - U k ). (3) 

The slower the velocity relaxes towards its desired value the smaller the amount of 
perturbations introduced in the system. On the other hand choosing a small value of A in- 
creases the number of iterations to reach equilibrium and thus decreases the computational 
efficiency Here we choose as a relaxation parameter A = 0.1. 



2.2 Lattice Boltzmann model 

The continuum hydrodynamics are described by the incompressible Navier-Stokes equa- 
tions 

— + (u-V)u = - Vp/p + uV 2 n + g, (4) 
V-u = (5) 

where u is the fluid velocity, p the pressure, p the density and g a body force. We use g 
to enforce Dirichlet boundary conditions, see discussion below. 

We solve the equations of motion (jlj) and © by using a Lattice Boltzmann algo- 
rithm [T3|. This approach follows the evolution of particle distribution functions fi on 
a d- dimensional regular lattice with z links at each lattice point r. The label i denotes 
velocity directions and runs between and z. DdQz + 1 is a standard lattice topology 
classification. The D3Q15 lattice we use here has the following velocity vectors v^: (0, 0, 0), 
(±1, ±1, ±1), (±1, 0, 0), (0, ±1, 0), (0, 0, ±1) in lattice units. 
The Lattice Boltzmann dynamics are given by 

/<(r + AtVi, t + At) = f^r, t) + - (/f (r, t) - f^r, t)) + At 9i (6) 

T 

where At is the time step of the simulation, r the relaxation time. The equilibrium 
distribution function f^ q is a function of the density p and the fluid velocity u defined as 

z z At 

P = J2fo ' P U = f^i + -^-8- ( 7 ) 

i=0 i=0 / 
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The equilibrium distribution function is chosen as 



1 " — ) w iP I — — + ,4 v * ■ S- ( 9 ) 



^ (r ' t)= ^( ^ — () 

where c s = l/v3 is the speed of sound and iUj are weights chosen as wo = 4/9, ti>j = 1/9 
for i = 1 — 6 and iUj = 1/72 for i = 7 — 14. The forcing term is defined as 

1 \ fvi - Ui (Vi ■ Uj) 

Performing a Chapman-Enskog expansion on the LB dynamics ^B] shows that equations 
(jlj) and (jHJ) are recovered with a kinematic viscosity expressed as 

„=(£>!l( T -i) (io) 

At 3 V 2 ; v ; 

where Ar is the lattice spacing. 

We enforce Dirichlet boundary conditions using a local forcing term g. The governing 
equation Equation (jlj) is rewritten as 

— = = -(u- V)u- Vp + vW 2 u = RHS (11) 

where u* is the velocity at time t + At with no forcing term considered. Including a forcing 
term leads to 

u'M + ^-u(r,^ MS + gM) (12) 

where u d is the desired velocity. Subtracting equations (fTTj) from (fT2~j) leads to an expression 
for the forcing term 

gM) = » J M + ^)-u'(r, t + A t ) (u) 

It is worth noting that evaluating u* within an LB method consists only of performing 
a streaming step. Indeed by construction of the equilibrium distribution function (pu = 
J2 fi 9 Vi, see [TO] for details), the second term of the right hand side of equation (jBj) does 
not change the velocity u. 



2.3 Hybrid model 

We use a domain decomposition algorithm to couple an MD description of a dense fluid 
with an LB model solving the incompressible NS equations. The computational domain is 
decomposed into two overlapping regions of an LB domain and an MD subdomain. Fig. Q 
shows the Schwarz decomposition used to converge to a steady state solution by alternating 
iterations between steady state solutions in the LB domain and MD subdomain. 
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Figure V. Domain decomposition. The dark gray area is computed within the MD sub- 
domain and (VGC) enforced within the LB domain or (VC) only enforced along a strip 
depicted by the dashed square. Converged solution is obtained by alternating iterations in 
the LB and MD domains. 

A Schwarz iteration t c consists of computing the continuum velocity field u c (t c ) with 
boundary conditions set by the previous atomistic cycle u a (i c — 1) and an external boundary 
condition that depends on the system considered. Then, u c (t c ) is used to set the boundary 
condition for computing u a (t c ). 

MD velocities are sampled in cells of same size as in the LB domain and are enforced on 
the continuum according to two coupling methods using equation (JT3J). The first approach 
corresponds to the one used by Werder et al. and is to impose MD velocities within a 
one cell wide strip located at a distance S s of the MD subdomain, see fig. ^ This velocity 
coupling (VC) approach does not enforce velocity gradients implying that the geometry 
of the system and the external boundary conditions dictate whether the hybrid solution 
evolve into a good approximation of the reference solution. To alleviate this issue we 
propose an enriched coupling method which enforces velocities, and implicitly velocity 
gradients (VGC) by imposing MD velocities on every common cell except within a strip 
of width 5 S close to the boundary (see fig. Q). We shall observe below that using the 
VGC approach leads to closer match with MD reference solutions than when using the VC 
method. 

We use the algorithm proposed by Werder et al. [13J to impose non-periodic boundary 
conditions (NPBC) on the MD system. Details of the algorithm can be found in [T3] . 

3 Results 

We first consider the flow past a CNT embedded normal to the flow direction in order to 
compare our results with We then discuss the flow through a short CNT embedded 
parallel to the flow direction. 
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3.1 Flow past a nanotube 

We apply the hybrid LB-MD algorithm to the case of the flow of argon around a CNT 
centered along the z-axis within a 30 x 30 x 4.254 nm 3 domain Q. The CNT is of chirality 
(16,0) with a radius of r = 0.625 nm. We choose the density of argon pA r = 1008 
kgm -3 and the temperature T = 215 K. This corresponds to the dimensionless state point 
(T*,p*) = (1.8,0.6) where T* = k B Te~ A ). Ar A and p* = pa 3 ArAr m A lA. We let k B be the 
Boltzmann constant, rriAr — 0.03994 kgmol -1 the atomic mass of argon, and A Avogadro's 
number. 

The MD sub domain size is 10 x 10 x 4. 254 nm 3 centered around the CNT which is 
subdivided into 20 x 20 x 1 sampling cells where 6465 argon atoms are initially equilibrated 
for 0.2 ns. The width of the strip around the boundary on which both MD and LB are 
simulated is 5 S = 2.5 nm. 

We consider an LB domain of size 60 x 60 x 1 covering the entire domain where lattice 
nodes are centered on the corresponding MD sampling cells. The viscosity of the LJ fluid 
is a parameter of the LB model and set to v = 0.745 • 10~ 7 m 2 s _1 We have performed 
a sensitivity analysis by increasing and decreasing the viscosity by 5% and found it to be a 
robust parameter with respect to accuracy. Dirichlet boundary conditions Mqo = u x = 100 
ms _1 are imposed at the inlet x = nm and outlet x = 30 nm. This high velocity is 
chosen to reduce the number of sampling iterations. The temperature is controlled by 
using a Berendsen thermostat ^H] with a time constant tt = 0.1 ps. We apply it cell-wise 
in all directions in the boundary cells where the velocity is prescribed and only in the 
z-direction in other cells. The hybrid model is run for 100 cycles which consists of running 
the LB simulation for 7 ns (15000 iterations) followed by an MD step equilibrating for 0.2 
ns (20000 iterations) and sampling for t s = 0.4 ns. We discuss below how t s affects the 
convergence of the results. 

Fig. |21 shows a comparison between hybrid converged solutions and an MD reference 
solution over the entire domain. The latter involves 58198 argon atoms and the temperature 
is controlled as in the MD subdomain. The reference MD velocity field is sampled over 20 
ns. A qualitative match with the reference MD solution is obtained when using the VC 
approach. The match between the hybrid solution and the reference MD solution becomes 
quantitative when using the VGC method. This is due to the fact that velocity gradients 
are not let free but implicitly imposed. A consequence of this is that discrepancies in the 
wake and on the sides of the CNT that appears when using the VC approach, and also 
reported in ^3], are not observed when using the VGC method. 

Fig. El shows the evolution of the norm of the velocity along the y — 15 nm and 
x = 15 nm plane. The system starts from a constant initial condition u x = Uoo = 100 
ms" 1 . This unphysical state leads to the formation of high velocity regions around the 
CNT. They then slowly spread out on the sides of the tube and increase the side velocity. 
Considering the VGC approach, we observe that the velocity upstream and downstream 
adjusts smoothly to the reference solution. The situation is different when employing the 
VG method. Indeed the hybrid velocity profiles show deviations from the MD reference 
solution especially downstream the CNT. It is worth noting that we get identical results 
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Figure 2: (left) Converged hybrid and (right) reference MD solutions of a flow past a 
CNT. (top) Velocities and (bottom) velocities and velocity gradients are enforced. The 
colorscale and contour lines depict the norm of the velocity expressed in ins -1 . Gray lines 
are streamlines. The thick squares represent (solid) the boundary of the MD subdomain 
and (dashed) the boundary of the overlapping region. 
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when initializing the system with initial conditions closer to the steady solution. 

In agreement with we observe that argon flow past a CNT normal to the flow 
results in a vanishing velocity at the surface of the CNT. In this particular configuration, 
the velocity boundary condition around the tube can therefore be approximated by a no- 
slip condition allowing us to solve the NS equations. Fig. El shows LB velocity profiles 
where the CNT is modeled off-lattice using an Immersed Boundary technique [20J. We 
will see however that this approximation is not applicable when we consider a CNT with 
a different orientation. 

We quantify the convergence towards the reference MD solution by defining an error e J 
between the hybrid solution at cycle j and the reference MD solution as 



e 



j 



where N is the number of cells in Q, and vlu,md are respectively the hybrid and reference 
MD velocities at cycle j in the cell k. 

Fig. E] shows the time evolution of the error which rapidly decays during the first 
10 cycles. The error then fluctuates around an average value which is a function of t s . 
Considering the VGC approach, we measure an average error between cycle 50 and 100 
of 1.3% for t s = 0.4 ns and t s = 0.8 ns and an average error of 1.9% for t s = 0.2 ns. We 
observe in fig. HJthat considering short t s leads to undesirable fluctuations whereas long t s 
decreases the computational efficiency of the model. An optimum is found to be t s = 0.4 
ns assumed hereafter if not otherwise specified. Fig. 0] also shows the time evolution error 
when using the VC method. The average error between cycle 50 and 100 is of 3.8% and is 
comparable to the the error evolution reported in |13j . 

The point-wise error between the hybrid and the reference MD solutions is depicted 
in fig. HJ Using the VGC technique leads to a localized and maximum 3.9% error region 
showing only on one side of the CNT. Everywhere else the error is up to 2.5%. The error 
is higher when using the VC approach where we observe a large error of 40% close to the 
CNT and up to 12% in the wake. 

The supercritical state point (T*,p*) = (1.8,0.6) is associated with less structural 
correlations than in a liquid state. In order to assess their effect on the convergence 
we have considered the liquid state point (131.fr, 1361 kgm~ 3 ) = (1.09,0.81) related to 
the kinematic viscosity v = 1.42 • 10~ 7 m 2 s _1 For similar parameters as for the 

state points (1.8,0.6), we measure an error of 2.0% between the hybrid solution and the 
associated reference MD solution. The error is slightly higher than the one obtained in 
the supercritical state and indicates that the LB-MD model can also be applied when 
considering liquid states. 



3.2 Flow through a nanotube 

We now consider the flow through a CNT driven driven by a constant velocity u x = 100 
ms _1 enforced at x = nm and x = 28. The CNT is centered along the x-axis in a domain 
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Figure 3: Evolution of the hybrid profiles of the norm of the velocity in the planes located 
(left) at x — 15 nm and (right) at y — 15 nm. (top) Velocities and (bottom) velocities and 
velocity gradients are enforced. The solid thick line is the reference MD solution whereas 
the solid line is the continuum solution where no-slip boundary condition is considered 
around the CNT. Hybrid solutions are plotted after various number of cycles. Dashed 
lines represent the boundaries of the MD subdomain. 
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Figure 4: (top) Evolution of the error e- 7 between the hybrid solution and the reference 
MD solution. Different sampling time t s and coupling techniques are considered, (bottom) 
Point-wise error according to equation (j!4|) shown in percentage and averaged between 
cycle 50 and 100. The sampling time t s = 0.4 ns. (left) The VGC and (right) the VC 
techniques are used. 
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Figure 5: Reference MD (top) and hybrid (bottom) solutions of the flow through a CNT 
driven by a constant velocity u x = 100 ms" 1 enforced at x = nm and x = 28 nm. 
The norm of the velocity is plotted. Each column corresponds to a different x — y plane 
at z = 8,4 and nm. Contour lines are plotted and expressed in ms -1 . Bold and thin 
rectangles depict the CNT and the MD subdomain respectively. 

of size 28 x 16 x 16 nm 3 and it is of chirality (16, 0) with a length I = 2.34 nm. We choose 
the dimensionless state point (T*,p*) = (1.8,0.6). 

The MD subdomain size is 14 x 10 x 10 nm 3 centered around the CNT and is subdivided 
into 14 x 10 x 10 sampling cells. 21278 argon atoms are initially equilibrated for 0.2 ns and 
5 S = 3 nm. The temperature is regulated by applying a Berendsen thermostat cell-wise in 
every direction within the boundary cells. We consider a 28 x 16 x 16 LB domain covering 
the entire domain where lattice nodes are centered in corresponding sampling cells. MD 
and LB domains are coupled via the VGC approach. The hybrid model is ran for 50 cycles 
consisting of running the LB simulation for 13.4 ns (15000 iterations) followed by an MD 
step equilibrated for 0.2 ns and sampled for t s = 0.4 ns. 

Fig. El shows a comparison between the hybrid and the reference MD solution. The 
latter consists of 108943 argon atoms sampled for 20 ns where thermal boundary conditions 
are as in the MD subdomain. There is an overall good quantitative agreement between 
the solutions. We observe regions above and under the CNT where the hybrid velocity 
is slightly faster than in the reference MD solution. The difference is at most 4%. We 
quantify the overall agreement by measuring average errors following equation ([T4*|) over 
different regions and during the last 30 cycles and obtain a global error e g = 2.6% over the 
whole domain and a local error e/ = 2.1% in the tube. 

The rather small diameter of the CNT leads to a higher argon density within the tube. 
We measure a 7% increase that can not be described by the incompressible NS equations. 
This issue could however be alleviated by considering a compressible LB model |2*T] . 

In this conformation, the velocities around the CNT are non-zero. The corresponding 



11 



continuum velocity boundary condition is unknown and a function of the system parame- 
ters PHI- A continuum solution approximating the reference solution is therefore unfeasible. 

The computational efficiency of the LB-MD model is estimated by comparing the time 
needed to compute one iteration of the reference and hybrid solution, respectively. We 
measure t re f = 0.56 s and th y b = 0.08 s in the case of the flow past a CNT. The hybrid 
solution is computed R = 0.56/0.08 = 7 times faster than the reference solution. We get 
a ratio R = 1.04/0.43 = 2.4 in the case of a flow through a CNT. Ratios R are of the 
same order as the volume ratio between the MD domain and subdomain. Small systems 
have been chosen here in order to compare hybrid to reference solutions. Considering 
larger systems would exhibit much higher ratios. For example, nanodevices (cx 10 3 nm 3 ) 
embedded in microscale systems (oc 1 /im 3 ) would lead to ratios of the order of 10 6 . 

4 Conclusion 

We have presented a hybrid model coupling a Lattice Boltzmann solution of the incom- 
pressible Navier-Stokes equations to a Molecular Dynamics Simulation of a dense fluid, 
using a Schwarz domain decomposition technique. The two descriptions are coupled via an 
exchange of velocities (VC) or velocity gradients (VGC). The applicability of the method 
was demonstrated in flows of liquid argon past Carbon Nanotubes normal and aligned with 
the flow velocity. 

We have first computed the flow past a CNT with an axis normal to the flow velocity. 
Using the VGC approach we observed quantitative agreement between the hybrid and the 
reference MD solutions with an average error of 1.3% that provides an improvement over 
previously reported results [T3] for the same configuration. We attribute this improvement 
to the fact that we enforce velocity within a region rather than along a strip implying 
that velocity gradients and in turn shear forces are implicitly exchanged between the two 
descriptions. We also show that the velocity vanishes around the tube in this conformation 
making a continuum solution a good approximation of the reference solution. This is not 
the case however when considering the flow through a CNT. The transport of argon in the 
CNT is captured within 2.1% with our LB-MD model. To the best of our knowledge, this 
is the first time that fully 3D hybrid simulations have been reported in the literature. 

In addition of being a tool to help investigating multi-scale physics by considerably 
reducing an otherwise prohibitive computational time, this approach largely extends the 
applicability of LB models by providing us with fully microscopic boundary conditions. 
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